Multiple Histogram Method for Quantum Monte Carlo 
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An extension to the multiple-histogram method (sometimes referred to as the Ferrenberg- 
Swendsen method) for use in quantum Monte Carlo simulations is presented. This method is shown 
to work well for the 2D repulsive Hubbard model, allowing measurements to be taken over a continu- 
ous region of parameters. The method also reduces the error bars over the range of parameter values 
due the overlapping of multiple histograms. A continuous sweep of parameters and reduced error 
bars allow one to make more difficult measurements, such as Maxwell constructions used to study 
phase separation. Possibilities also exist for this method to be used for other quantum systems. 

PACS Numbers: 02.70.Lq,71.10.Fd,02.50.Ng,02.60.Gf 



I. INTRODUCTION 



II. SINGLE HISTOGRAM METHOD 
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When making calculations using the Monte Carlo 
method, one often would like to make measurements of 
some observable as a function of the parameters of the 
hamiltonian. To do this the standard procedure is to 
perform a run at one setting of the parameters until a 
measurement of the observable with sufficiently small er- 
ror bars is produced. One then moves on to another 
setting of the parameters, and so on until a large discrete 
set of measurements of the observable is produced. This 
can require large amounts of computer time depending 
on the details of the hamiltonian and the Monte Carlo 
method one has chosen. Building an approximation of a 
continuous function for the desired observable point by 
point would require an even larger set of discrete points, 
and a related increase in computer time. 

Continuous functions of observables are useful for in- 
numerable activities, such as verifying functional depen- 
dencies and looking for phase separation using-.Maxwell 
constructions .El The multiple histogram methodj (MHM) 
allows one to produce these continuous functions for a 
classical hamiltonian. In this paper, it will be shown 
that under certain circumstances the MHM can also be 
applied to quantum mechanical hamiltonians as well, for 
example the two-dimensional Hubbard model.™ Observ- 
ables of this hamiltonian will be measured-,using standard 
quantum Monte Carlo (QMC) techniques^ and the MHM 
will be applied to obtain observables as continuous func- 
tions of the parameters of the system. As in the classical 
case, the use of overlapping histograms reduces the errors 
below those of a single measurement. Furthermore, the 
MHM will be applied specifically to the density and en- 
ergy of the Hubbard model as functions of the chemical 
potential to search for signs of phase separation. 



In order to obtain a variable as a continuous function 
of thermodynamic variables like the inverse temperature, 
ft = (k^T)^ 1 , and the chemical potential, /i, one can use 
a technique originally, used by Salsburg et al.p extended 
by Valleau and (fj*fd,Q and further extended by Swendsen 
and Ferrenberg.ErQ To see how this can be done for the 
chemical potential, let us assume we have a classical 
hamiltonian, 



Ti. = E — pn, 



(1) 



where E is the energy and n is the particle density. From 
this hamiltonian, one can derive a partition function for 
a particular value of p, 



Z^ = Y^ P( n ) e 



(2) 



where p(n) is a density of states which for these pur- 
poses we can leave undetermined. We can also derive a 
probability to have a given density for this value of 



p(n)e 



(3 fj,n 



(3) 



Substituting Eqn. ||into Eqn. ^ allows us to find the prob- 
ability for a value of n at a different value of the chemical 
potential, \i! , 



PA") 



P (1 (n)e«^' 1 '" 



(4) 



Eqn. | can be extended for use with many histograms 
taken at varying parameters, as discussed by Ferrenberg 
and Swendscn.o Using multiple histograms allows one to 
increase the range of validity of the technique. 
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FIG. 1. A histogram of the particle density sampled af- 
ter each sweep in a Quantum Monte Carlo simulation of the 
Hubbard model on a 8 by 8 lattice with fit = 8, U/t = 8, and 
fj, = —1.10 and —1.30. One can see well defined peaks at par- 
ticle densities corresponding with integer numbers of holes. 
It is also possible to see how the two histograms at different 
parameters overlap each other, allowing one to be normalized 
with respect to the other. 

To demonstrate the use of this method for quantum 
hamiltonians, the Hubbard hamiltonianB will be used 

iv) <?={U} 
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where t is a parameter to set the strength of electron hop- 
ping (kinetic energy), c^J is the annihilation (creation) 
operator for an electron at site % with spin a = {Tl}, U 
is a parameter to set the strength of on-site coulomb re- 
pulsion, /i is the chemical potential, and hi a — c icr Ci a . 
Clearly there are quantum operators in the Hubbard 
hamiltonian, so the classical derivation given earlier for 
the probability as a function of fi will only hold if the 
system happens to be in states of definite particle num- 
ber. In determinantal QMC, a Hubbard-Stratonovich 
(HS) transformations leads to a bilinear fermion form 
for the interaction. One sums over all possible fermion 
states and then over the HS field. Thus a given HS con- 
figuration comes from a trace over all particle numbers 
so that it is not guaranteed that a given HS configura- 
tion will have a definite particle number. It is this fea- 
ture that distinguishes this problem from the classical 
statistical mechanics problem. However if one finds that 
the various HS configurations are indeed characterized by 



an integer fermion occupation, then as we will discuss, 
one can proceed. Here for the 2D Hubbard model near 
half-filling, the charge gap provides an energy barrier to 
non-integer filling. As can be seen in Fig. |lj, when one 
runs simulations at a sufficiently low temperature, one 
observes peaks in the histogram corresponding to states 
with integer numbers of holes. If these states truly have 
a definite number of particles, the peaks should scale us- 
ing the grand canonical distribution for runs at different 
values of /x. In Fig. ||, such a test is performed and the 
peaks do scale with the grand canonical distribution as 
desired. 
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FIG. 2. Normalizing the peaks in histograms such as the 
ones shown in Fig. jj], yield the probability for a certain parti- 
cle density. One can then change the fj, used in the simulation 
and follow the relative probability of this peak as a function 
of fj,. If the states truly have a definite particle density they 
should follow a curve given by the grand canonical distribu- 
tion which is shown in the solid line. 

Note that in Fig. [|, multiple histograms are presented 
in order to emphasize the fact that the MHM allows one 
to obtain even more information when several histograms 
are used. By normalizing each histogram with respect to 
the others, a continuous function of the desired observ- 
able may be obtained. In the results that follow, we have 
used many overlapping histograms in order to fully cover 
the range of simulation .variables of interest. 

In parameter regimeso where these peaks are present 
in the histograms of particle number, one can proceed to 
take any operator that conserves particle density and de- 
termine its behavior as a continuous function of \i. This 
can be done by using a reduced number of simulations at 
different values of fi, as long as the respective histograms 
of particle density overlap. To prevent unwanted overlap 
from very distant histograms, it was found necessary to 
introduce a Gaussian cutoff. The cutoff reduces the ef- 
fective weight of a histogram as the difference between 
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the value of \i for that histogram and the value of // at 
which the MHM method is run increases. As an added 
benefit, the error bars for this continuous function of /x 
will be reduced from those obtained from any one simula- 
tion, because one can use information about the operator 
from all of the simulations performed. 




0.92 

-1.85 -1.6 -1.35 -1.1 -0.85 -0.6 



FIG. 3. The original points and the MHM method continu- 
ous curves are shown for the particle density (n) vs. chemical 
potential (fi) for two values of the inverse temperature, (3. 
All simulations are performed on an 8 by 8 square lattice at 
U/t = 8. Error bars on the original data points are shown as 
bars with hats originating from the points, and error bars on 
the MHM line are shown as bars without hats originating at 
regular intervals from the dashed line. 

III. RESULTS 

In Fig. |^ one can see what happens when the Multiple 
Histogram Method (MHM) is applied to a set of his- 
tograms measuring the particle number during a QMC 
simulation of the Hubbard model. The chemical poten- 
tial, /x, is varied over a range starting from half-filling 
(zx = and thus n — (n) — 1) through increasingly nega- 
tive values of /x. The individual runs are shown as points 
with associated error bars, and the MHM is applied to 
get a continuous function of n vs. fi with error bars 
shown as vertical bars around the dashed line. The error, 
bars are checked through both the bootstrap methodtHl 
and by incorporating the errors due to correlations in the 
reweighted samples and their finite size.tiJ To see how the 
line changes with different values of the inverse temper- 
ature, f3, the process is repeated at fit = 4. 

As one can see, the error bars on the MHM line are 
comparable to those from the original points, and the 
MHM continuously fits the original points. As the simii=. 
lation temperature decreases (large fi) the sign problems 
becomes much more severe. For the histogram method 
this means that while there are peaks in the histograms 



for n vs. fi for configurations with positive signs, they are 
wiped out by nearly identical peaks with negative signs. 

In Fig. ||, the MHM is applied to another interest- 
ing observable, the equal time antiferromagnetic spin 
correlation function, S zz (ir, 7t), the Fourier transform of 
((riij — riii)(njj — "ji)} taken at momentum (-7T, tt). As 
the system moves closer to half-filling (/i =^-,0), the anti- 
ferromagnetic order increases as expected.Ej 



♦ Orig pt=4.00 
MHM pt=4.00 
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FIG. 4. Original points and the resulting MHM method 
continuous curve are shown for the equal time zz antiferro- 
magnetic spin correlation function (S zz (tt, tt)). As the pa- 
rameters are moved closer to half-filling (/i = 0), the system 
begins to develop antiferromagnetic order as shown by the 
increasing value of S zz (n,n). The simulations are performed 
on an 8 by 8 square lattice at U/t — 8. Error bars on the 
original data points are shown as bars with hats originating 
from the points, and error bars on the MHM line are shown 
as bars without hats originating at regular intervals from the 
dashed line. 



IV. PHASE SEPARATION 

Recently, computational evidence for the existence of 
phase separation in the strongly correlated t-J and Hub- 
bard models, bath for and against, have been presented in 
the literature.ll-a One of the best ways to search for phase 
separation is through the use of a Maxwell construction.!!! 
A Maxwell construction consists in its simplest form as 
a graph of energy vs. density. Thermodynamics wants a 
system to minimize its free energy, and so a phase sepa- 
rated region will appear on the graph as a region of up- 
wards curvature (§-4r < 0). In such a region the energy 
of the phase separated state will be lower than that of the 
homogeneous state, and phase separation is thus favored. 
It is of interest, therefore, to apply the MHM to obtain 
continuous energy versus density measurements in order 
to investigate the possible existence of phase separation. 
In Fig. [|, a Maxwell construction for the 2D Hubbard 
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model at j3t = 4 is shown. With just the original data 
points it is nearly impossible to say anything definitive 
about the curvature of the line, however when the MHM 
method is employed, the lack of any clear signs of up- 
wards curvature becomes much more apparent. Thus it 
is possible to conclude that one does not see phase sep- 
aration for this range of parameters, a conclusion that 
would have been difficult to make without using the full 
information from the histograms. 




-0.99 1 1 1 1 d 

0.96 0.98 1 

n 

FIG. 5. Original points and the resulting MHM method 
continuous curve are shown for a Maxwell construction of en- 
ergy vs. particle density (n); some points have been omitted 
for clarity. Phase separation would be indicated by a region 
of the curve with upwards curvature, which appears in nei- 
ther curve. Error bars for both the uncertainty in energy and 
particle density on the original data points are shown as bars 
with hats originating from the points, and error bars on the 
MHM line are shown as dotted bars without hats originating 
from the dashed line. 



V. CONCLUSIONS 

In some circumstances, the MHM method can be used 
effectively for hamiltonians containing quantum opera- 
tors, such as the Hubbard hamiltonian just described. 
The ability to generate continuous functions of an observ- 
able makes a number of analysis techniques much more 
feasible. For_ example, in order to perform a Maxwell 
construction,!! a continuous function of energy as a func- 
tion of particle density is desired. Previously a func- 
tion like this would have been constructed laboriously 
one point at a time using the averages of large simu- 
lations, and conclusions drawn using reasonable guesses 
about what function would fit the calculated points. Us- 
ing the multiple histogram method, the large quantity 
of information present in the histograms of each simu- 
lation is brought to bear using simple statistical physics 



arguments. The information from the simulations is thus 
used more efficiently and smaller error bars and reason- 
able continuous functions are produced. 
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